Lyman-a transfer in primordial hydrogen recombination 
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Cosmological constraints from the cosmic microwave background (CMB) anisotropics rely on 
accurate theoretical calculations of the cosmic recombination history. Recent work has emphasized 
the importance of radiative transfer calculations due to the high optical depth in the H I Lyman 
lines. Transfer in the Lya line is dominated by true emission and absorption, Hubble expansion, 
and resonant scattering. Resonant scattering causes photons to diffuse in frequency due to random 
kicks from the thermal velocities of hydrogen atoms, and also to drift toward lower frequencies 
due to energy loss via atomic recoil. Past analyses of Lya transfer during the recombination era 
have either considered a subset of these processes, ignored time dependence, or incorrectly assumed 
identical emission and absorption profiles. We present here a fully time-dependent radiative transfer 
calculation of the Lya line including all of these processes, and compare it to previous results that 
ignored the resonant scattering. We find a faster recombination due to recoil enhancement of the 
Lya escape rate, leading to a reduction in the free electron density of 0.45% at z = 900. This results 
in an increase in the small-scale CMB power spectrum that is negligible for the current data but 
will be a 0.9o" correction for Planck. We discuss the reasons why we find a smaller correction than 
some other recent computations. 

PACS numbers: 98.70.Vc, 98.62.Ra, 32.80.Rm 



I. INTRODUCTION 

The cosmic microwave background (CMB) anisotropy 
has proven to be one of the most useful and robust cos- 
mological probes. The most recent example has been the 
Wilkinson Microwave Anisotropy Probe (WMAP) satel- 
lite, which has provided key information on the com- 
position of the high-redshift Universe, the distance to 
the smface of last scattering (via the acoustic peak po- 
sition), cosmic reionization, and the primordial power 
spectrum P, i, The development of precision cos- 
mology with the CMB will continue with the upcoming 
launch of the Planck satellite. Planck data will be essen- 
tial for dark energy studies: its measurement of ilmh^ 
and Da{z = 1100) will be needed to break degeneracies 
with 'w{z) in supernova data and set the length of the 
baryon acoustic oscillation standard ruler, and the pri- 
mordial power spectrum will be needed to interpret tests 
of the growth of structure using weak lensing and clus- 
ters 0, S @1 • The measurement of the primordial spectral 
index rig will be essential for constraining models of in- 
flation; indeed it is the precision measurement of rig from 
WMAP (and upper limits on tensors )_that have opened 
the era of tests of inflaton potentials [1, 0] • 

In addition to its versatility, one of the advantages of 
the CMB as a cosmological probe is that the generation 
of CMB anisotropics is understood from first principles: 
they arise due to linear perturbations on a homogeneous. 
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isotropic background of baryons, photons, dark matter, 
and neutrinos. There are now several public codes that 
solve these equations and compute the angular power 
spectrum Ci and whose numerical accuracy has now 
been checked to 1 part in 10^ [l3|- However the photons 
and baryons interact via Thomson scattering, so these 
codes require as input the ionization fraction Xe{z). This 
is a complicated non-equilibrium physics problem as one 
must track the level populations of recombining hydro- 
gen and helium atoms, as well as following the effects of 
the emitted radiation and its re-absorption. 

The first attempts to sol ve p rimordial recombination 
were carried out by Peebles 11 1 and Zel'dovich et al [l^l , 
and the progress of CMB experiments has driven theo- 
rists to construct ever more accurate models of recombi- 
nation. Most CMB codes in common use today obtain 
Xeiz) from the Recfast code by Seager et al [H, 
which is based on the multi-level atom (MLA) approxi- 
mation. The MLA follows the level populations of each 
type of atom considered (H i. He i. He ii) under the in- 
fluence of bound-bound and bound-free transitions. The 
optically thick lines (the H i Lyman series and its He i 
and He ii analogues) are treated using the Sobolev ap- 
proximation [lij to avoid the need to track the radiation 
field explicitly. Two-photon decays from the H i 2s, He i 
2^5*0, and He ii 2s levels and their inverse process (two- 
photon absorption) are important and also included. 

The picture that emerged from this work is one in 
which hydrogen recombination occurs far more slowly 
than the Saha equation predicts. The reason is the pro- 
duction of radiation: a hydrogen atom that recombines 
directly to the ground level (Is) emits a Lyman contin- 
uum photon. Once ^ 1 out of every 10^ H atoms has 
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recombincd the Universe becomes optically thick to this 
radiation and each direct recombination is immediately 
followed by an ionization. Recombination can therefore 
occur only indirectly, via recombination to the H i ex- 
cited levels {nl, n > 2) and radiative cascade to Is. But 
here another bottleneck occurs: the buildup of radiation 
in the Lya {2p Is, 1216A) line. In the absence of a 
sink for Lya photons, a hydrogen atom in the 2p level 
can only reach the ground level by exciting another hy- 
drogen atom to 2p. The intense thermal CMB radiation 
during the recombination epoch is easily capable of ion- 
izing a hydrogen atom from the n = 2 level, so again the 
process produces no net recombinations. The two major 
processes that break the bottleneck are (i) the cosmologi- 
cal redshifting of photons out of the Lya line, and (ii) the 
2s ^ Is two-photon decay (the photons have a contin- 
uous energy distribution and hence are not immediately 
re- absorbed in hydrogen lines). Similar physics applies at 
earlier epochs to helium recombination, albeit with some 
additional complications due to the more complex energy 
level structure. 

In the past several years, many authors have re- 
investigated the recombination problem with attention 
to small physical effects, in particular neglected pathways 
that might allow hydrogen atoms to reach the Is level. In 
particular, the Planck mission will require much better 
accuracy than was desired in 1999 when the original ver- 
sion of Recfast was written. The new effects considered 
have included two-photon decays from the ns and nd lev- 
els of H I JT^JTzlJlsl, [H, iOl and the n^S and n^D levels 
of He I |l6l. ITI. l2l|: H i continuum absorption of He i 
584A line radiation [l^, [H, ; stimulated two-photon 
decays and two-photon absorption [l^, [2^, Raman 
scattering [l^, [2l| ; resolution of the H i ^-sublevels [l^l ; 
and forbidden transitions in He i [l3, [H, [l^, [l^ . In some 
cases the corrections to the CiS were > 1%, sufficient to 
bias the Us measurement from Planck by several sigma 
[29| . Some of these corrections are now incorporated into 
the most recent version of Recfast [s^I • 

An additional correction is the deviation from strict 
Sobolcv behavior in the very optically thick H i Lya 
line {2p Is). This is the subject of this paper. The 
Sobolcv approximation is based on the assumptions of (i) 
identical absorption and emission profiles, (ii) complete 
frequency redistribution in each line scattering, (iii) ab- 
sence of any other absorption or emission processes active 
in the same frequency range as the line, and (iv) quasi- 
stationarity, i.e. that the level populations and radiation 
field change little during the time it takes for a photon to 
rcdshift through the line. None of these approximations 
arc quite valid during the recombination era. Previous 
authors have relaxed some of these assumptions individ- 
ually, but there exists no comprehensive treatment. Kro- 
lik ll^l considered the partial redistribution during 
Lyman-a scattering, and found only a small correction 
due to atomic recoil. Rybicki & dell'Antonio [s^ relaxed 
the quasi-stationary assumption and found that the re- 
laxation timcscalc for the Lya line was a factor of a few 



shorter than the recombination time; in the context of 
modern experiments this would probably imply a signif- 
icant correction to the recombination history. Hirata & 
Switzer [2lj argued in the context of the He i 584A line 
that the deviation of absorption versus emission profiles 
would lead to a somewhat increased escape probability, 
an effect which was verified to be important for H i Lya 
by Hirata [l^. There are also recent treatments of dif- 
fusion and atomic recoil by Grachev & Dubrovich [s^l 
and quasi-stationarity by Chluba & Sunyaev [s^ . How- 
ever there is not yet a treatment relaxing all of (i-iv) 
simultaneously. 

In this paper, we will first review the physics of the 
recombination era with an emphasis on the role of Lya 
escape (i jll)) . Next we describe a numerical method for 
solving the Fokker-Planck equation in the vicinity of Lya 
and grafting it on to an existing MLA code, and present 
the results obtained by this technique f ijIII[) . We then 
describe an analytic approach to the Lya escape, which 
is much simpler than the fully numerical technique and 
captures the essential physics f illVp . We describe impli- 
cations for precision CMB observations in ^|Vl We con- 
clude in WII Appendix [A1 describes computation of two 
special functions xi^i and T{W, S) introduced in this 
paper. 

For ease of comparison, we use the same fiducial cos- 
mology as in Ref. [H: n,nh'^ = 0.13, fJb/i^ = o.022, 
Tqmb = 2.728 K, helium mass fraction Y = 0.24, and an 
effective number of masslcss neutrinos N^, = 3.04. We 
draw heavily on the work of Hirata [l^ , where our MLA 
code was first presented. 



II. BACKGROUND 

In this section, we begin with a review of the basic 
definitions f ^II A[) . We then review the MLA method and 
its extension to two-photon transitions ( W B[) . We then 
describe the physical formulation of the Lya diffusion 
problem ( giIC|) . 



A. Definitions 

Since this paper extends the code of Hirata [l^, we 
begin with the same notation and basic equations as in 
that paper, and add new variables as needed for the 
Lya problem. The total density of hydrogen nuclei is 
flu oc where a is the scale factor. Ionization fractions 
are given by the free electron abundance Xe = ne/nn 
and free proton abundance Xp = np/nn relative to hy- 
drogen. Since this paper concerns the epoch after com- 
pletion of helium recombination and before any hydro- 
gen converts to molecular or intermediate (H~, HJ, etc.) 
forms, we have Xg = Xp. Photons are described by 
the phase space density f{E). The abundance of spe- 
cific energy levels of the hydrogen atom will be given by 
Xni = n\i{ i(nZ)]/nH. Level degeneracies = 2(2/ -I- 1), 
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energies Eni = —hlZ/n^, and Einstein coefficients Ay- 
have their usual meaning; we denote the hydrogenic Ry- 
dberg (in frequency units) by 7?., and set Boltzmann's 
constant equal to 1 so that temperature has units of en- 
ergy. 



B. Multi-level atom and tviro-photon transitions 

The Hirata [l^ code follows bound-bound and bound- 
free transitions involving each level of H i. The treatment 
of bound-free transitions and of the matter temperature 
is not altered by this paper and the corresponding equa- 
tions will not be repeated. For the bound-bound case, 
the rate equation is 



Xi\hh 



(1 + f]i+)X] - —fji+Xi 
9i 



■ Aij Pij 

j<i 
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-/y + Xj - (1 + f^J + )x^ 



(1) 



where the sums are over levels j that are above [j > i) or 
below (j < i) the energy of level i, and fji+ is the phase 
space density on the blue side of the line connecting levels 
i and j. The escape probability Pji is expressed in terms 
of the Sobolev optical depth Tji via 



P3^ = 



1 - e-^^' 



and the optical depth is 



9i' 



(2) 



(3) 



Of particular importance to us is the optical depth in the 
Lya (2p — * Is) line, which is typically of order 10^-10^. 

The phase space density fjij^ on the blue side of each 
line is required since these photons redshift into the line 
and can be absorbed or cause stimulated emission. In 
most cases it can be treated as a blackbody, but in the 
case of the Lyman series lines one must consider the non- 
thermal nature of the radiation field. In particular, pho- 
tons emitted in the Ly/3 line (1026A) will eventually red- 
shift into Lyg and begin exciting atoms. This "feedback" 
process 22, is taken into account by looking up the 
phase space density that emerged from Ly/3 at a previous 
epoch and taking this as input for the Lya calculation. 

Two-photon decays of the form 



H(nO-^H(ls)+7 + 7 



(4) 



produce photons at a rate per unit frequency per H atom 
of 



A„;(l/) 



dv 



(1 + fu){l + fu')Xnl - —fufv'Xls 
9ls 



where dA„; / dv is the spontaneous two-photon decay rate. 
One can easily incorporate the total decay rate 



Xnl\2i = - I /^rd{v)dv (6) 

/(l-)i-2)TC/2 

in the system of ODEs for the excited levels. However 
one must also take account of the radiation produced: 
if one of the emitted photons has an energy exceeding 
Lya, it will eventually redshift into the Lya resonance 
and excite a hydrogen atom. This is done by the vir- 
tual level method [l^, which implements Eq. ^ by in- 
troducing purely as a mathematical device a new virtual 
level of the hydrogen atom with energy Eis -f hv (where 
v is the frequency of the higher-energy photon) and in- 
finitesimal degeneracy. The decay in Eq. (j4]) can then 
be treated as a sequence of one-photon decays. The 
required choice of effective one-photon rate coefficients 
required for this procedure to work are given in §IVA 
of Ref. [l^; these choices also account for two-photon 
absorption. The same procedure also applies to other 
two-photon processes such as Raman scattering and two- 
photon recombination/ionization. Most importantly, the 
same feedback machinery used for one-photon transitions 
is automatically applied to the two-photon transitions. 

The two-photon differential decay rate dKni/dv con- 
tains resonances associated with the allowed 1-photon 
decays. For example, the two-photon decay rate from 
3d — !■ Is has a resonance corresponding to the sequence 

H(3d) ^ H(2p)-h7(Ha), H(2p) ^ H(ls)+7(Lya). (7) 

These resonances are physical since such decays actually 
do produce two photons [l^ • They do however result in 
very large decay rates if one naively computes the quan- 
tity 



A tot 
3d 



1 /-sK/grfA; 
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dv 



dv = 6.5 X 10^ 



(8) 



(5) 



(compare to 8.2 s^^ for the 2s level). Naively the exis- 
tence of such a rapid — > Is two-photon decay process 
would dramatically accelerate recombination, but since 
most of these decays produce photons within the Lya line 
this turns out not to be the case: almost every Sd Is 
two-photon decay is immediately undone by Lya absorp- 
tion. Ultimately the inclusion of these two-photon decays 
produces only 0(1%) corrections to the recombination 
history when one tracks the radiation field as well p^ . 
It is important to note that this result can only be ob- 
tained by consideration of radiative transfer in the Lya 
line 

Overall, this yields a system of ordinary differen- 
tial equations (ODEs) for the atomic level populations. 
These have explicit dependence on their history so the 
results of the integration must be stored for future refer- 
ence. Time steps are equally spaced in In a, with a fidu- 
cial choice of A In a = 4.25 x 10~^. The excited atomic 
levels are treated using the steady-state approximation. 



i.e. assuming that the atom reaches the ground state or 
is ionized in a time short compared to the recombination 
timescale. 



C. The Lyof diffusion problem 

The treatment of radiative transfer in the Lyman series 
described in ijllBI correctly describes many processes. It 
contains correct emission and absorption profiles for all 
of the Lyman lines, including interference between neigh- 
boring resonances, and is fully time-dependent. How- 
ever, one key piece of physics is missing: it neglects the 
Doppler shift due to the motion of the atoms. For con- 
tinuum processes far from resonance this is a minor er- 
ror. It is a major omission in the vicinity of the Lyman 
lines, where repeated scattering of photons off of hydro- 
gen atoms can cause them to undergo a random walk in 
frequency (frequency diffusion). This is especially true 
for Lya because of its very high scattering-to-absorption 
ratio. This section is devoted to a qualitative discussion 
of the physical processes involved; quantitative compu- 
tations arc deferred to mill 

The behavior of the radiation intensity near the Lya 
line is a competition among several effects. The resonant 
scattering of photons off of hydrogen atoms, 



virtual 



(9) 



allows photons to exchange energy with the kinetic de- 
grees of freedom of the matter but does not change 
the number of photons in the Lya resonance. There- 
fore if only this process were active it would drive the 
radiation intensity toward a modified blackbody fi, ~ 
|-g(/ii>-^i)/Tm -where the chemical potential is a con- 

stant. Since near Lya / ^ 1 the Bose-Einstein statistics 
of the photon are negligible and we expect oc e"'*''/^™. 
A second process is the Hubble expansion, which moves 
photons to lower frequency at a constant rate v = —Hv. 
A third process is the "true" emission and absorption 
of the Lya photons - that is, emissions and absorptions 
that are not part of the scattering process (Eq. [5]). In 
this case, the line profile approaches the form shown in 
Figure [TJ near line center the emission, absorption, and 
scattering are dominant. The emission and absorption 
set the phase space density of photons at line center to 
be that of equilibrium, /{vhya) — X2p/ (3xis), and the 
scattering then gives the frequency dependence 



3x 



Is 



(10) 



At some frequency z/trans sufficiently far to the red side 
of the Lya line, a transition occurs where the time to 
redshift out of the line becomes shorter than the time 
to diffuse back to line center. Beyond this point, the 
phase space density approaches a constant fhya--, which 
is related to the net flux of photons to the red side of the 



/v 



Lyu- 

__ 



3x, 




'Lya 



FIG. 1: A schematic representation of the photon phase space 
density near Lya. The dashed line shows the chemical equilib- 
rium solution, Eq. (|10|l . The true solution (solid line) comes 
to equilibrium near line center. As one moves to the red side 
of the line, the rate of absorption, emission, and scattering 
decrease, and at some point {v ~ Vt-cuns) the photons fall out 
of chemical equilibrium and redshift out of the line. 



Lya line via 



d(# photons) 87riJ 



dVdt 



X3 



fhyc 



(11) 



(Here SirH/Xf^^^ is simply the number of photon modes 
per unit volume per unit time that redshift through the 
frequency I'Lya — £•) Qualitatively, the rate at which 
photons redshift out of Lya is determined by Eq. PH)) at 
the transition frequency, combined with Eq. (jlip . This in 
turn gives the net 2p ^ Is decay rate (after a correction 
involving fhya+ for photons that redshift in to Lya and 
excite atoms is applied). In equations, we have 



8ttH 



X2p -/t(i/t„„,-i.Lvo)/rn, 



3a; 



Is 



hy 



a+ 



(12) 

This equation with I'trans ~ vj^ya and fhya+ = 
(•g'ji'Lyc./Tl. _ 2)^1 gives the Peebles [ll[ transition rate, 
equivalent to the Sobolev rate in the limit of r ^ 1; an 
excellent description of the physics can be found in §6 
of the book by Peebles [131 ■ The same equation, with 
various estimates of the "effective" t'trans, is the underly- 
ing conceptual reason for the accelerated recombination 
found by authors who considered repeated Lya scattering 
[33 |. (Recombination is accelerated rather than deceler- 
ated because f trans < J^Lya-) 

The role of emission and absorption in the Lya damp- 
ing wings must also be considered. It was at first ar- 
gued that because the emission and absorption profiles 
are the same (they are both Voigt profiles), this process 
tends to smooth out the frequency dependence of and 
suppress the boost factor e~''('^*'='°=~'^LyQ)/T„ -^^ ([T^ 

[32| . However, true absorption of Lya photons is actually 
a multiple-photon process: the virtual H(2p) atom must 
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not decay back to the ground state, but rather absorb 
another photon, e.g. 

H(ls) + 7(Lya) + 7(Ha) H(3s, M). (13) 

It follows that the Lya absorption profile actually de- 
pends on the color temperature of the ambient CMB 
near Ha, because of energy conservation: if the first pho- 
ton is from the red tail of Lya then the second must be 
from the blue tail of Ha and vice versa. In contrast, the 
emission profile, i.e. the reverse of Eq. (fT3|) . is simply 
the Voigt profile (with small corrections due to neigh- 
boring resonances and stimulated emission) . Thus in the 
cosmological context, the absorption profile is enhanced 
relative to the emission profile in the blue wing, and sup- 
pressed in the red wing, so true absorption and emission 
also tend to establish a red tilt to the radiation spec- 
trum. We can also understand this effect thermodynam- 
ically: since true absorption of a Lya photon followed by 
true emission (a process termed "incoherent scattering" 
by Krolik) exchanges the Lya photon's energy with the 
low-energy Ha photons, it follows that this process also 
tends to drive the radiation spectrum toward a modified 
blackbody, but this time with temperature T,-. At high 
rcdshifts we have Tr w Tm and so once again Eq. (fTO]) 
should apply if true emission and absorption are domi- 
nant processes. Their consideration leads to a spectrum 
similar to that of Fig. [T] and hence to a transition fre- 
quency and enhanced redshifting rate just as does Lya 
scattering. This effect has been seen in the consideration 
of two-photon decays; see e.g. §VB of Ref. [Toj . 

A fourth process is the 2s ^ Is two-photon contin- 
uum: a small fraction of this continuum overlaps the red 
damping tail of Lya (especially when modifications due 
to stimulated emission are considered) and this should 
be taken into account. 

A complete theory of recombination must take full ac- 
count of all these effects, and yield predictions for the 
2p — > Is decay rate and Lya radiation profile. It should 
go well beyond the conceptual discussion of "ftrans" and 
account for non-time-steady effects. 

The basic strategy in both our numerical and analyti- 
cal methods is to excise the region immediately surround- 
ing the Lya line from the aforementioned two-photon ra- 
diative transfer code and replace it with a method that 
solves the diffusion equation at high resolution. A hy- 
brid method is necessary because the resolution required 
near Lya (we must resolve the Doppler width of the line) 
would be too expensive if applied to the entire spectrum. 

III. NUMERICAL METHOD 

In our numerical method, we treat the photons near 
the Lya using a time-dependent Fokker-Planck method. 
This replaces the complicated redistribution of photons 
with a partial differential equation (PDE) that can be 
solved numerically. Wc first set up this equation and its 
method of solution, and then comment on the initial and 



boundary conditions and interface to the MLA code. The 
Fokker-Planck and MLA codes are interdependent, since 
the Fokker-Planck code requires boundary conditions and 
level populations from the MLA code, whereas the MLA 
code needs the Fokker-Planck code to determine the Lya 
decay rate and the flux of photons emerging from the 
red wing of Lya. We solve this problem by alternately 
running each code until convergence is achieved. An in- 
terface script passes key variables between the two codes, 
such as the phase space density of photons that redshift 
across the frequency boundaries. 

We divide the region of frequency near the Lya into 
bins spaced equally in In usually with spacing A In j/ = 
8.5 X 10~^. This is fine enough that it takes many steps 
for a photon to redshift through a Doppler width (the 
condition is that Aln;^ <C ^Tm/^iHC^)- We denote by 
Ni the number of photons per H nucleus per frequency 
bin in the i*^ bin, and the frequency associated with the 
i^^ bin by Vi. This is related to the phase space density 
via 

N. = ^4^/(^0. (14) 

This evolves with time according to Hubble redshifting, 
emission, absorption, and scattering: 

N^=N^\■a+N^\crn + N^U+N,Uc. (15) 

We use M = 2001 bins in our standard case, and place 
the lo'^ bin [where = |(M - 1) = 1000] at Lya line 
center. Then 

= ^'Lyaexp[(i - io)Alnz/]. (16) 

A. Transport in the Lya line 

Here we describe the computation of each term in 
Eq. nil). 

1. Emission 

We next consider the true emission of Lya photons. 
In the case of emission, the photons arise via two-photon 
decays from the n > 3 states, Eq. ([7]). In the limit where 
we are close enough to Lya line center to neglect the 
other resonances and the variation in photon phase space 
factors across the line, we may write this as 

= Xl^"r„^2p^'2p^i'''^v('^»)Ai'«. (17) 

u 

This is a sum over all higher-energy states of hydrogen u, 
where Tu~t2p is the rate per second at which those states 
decay to the 2p state, (/)v(i^) is the Voigt distribution, 
^2p^is is the branching fraction for 2p to decay to the 
ground level (P2p^is = 1 in vacuum), and Ah'i is the bin 
width. 
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Rather than keep track of all the variables necessary 
to make that calculation, we retrieve the Lya production 
rate in photons per hydrogen atom per Hubble time from 
the interface ( mil B[) . Equation (fTT]) simply reduces to 



NA 



Hn(|>Y{v^)Avi 



where 



(18) 



(19) 



denotes the Lyman a production rate. 

In reality the neglect of photon phase space factors and 
other resonances in Eq. (fT8)l is not correct to the desired 
accuracy. We have thus implemented a "corrected" form 
of the equation. 



NA 



HIl£{u,)^v{v,)Av„ 



(20) 



where Si^v) is the correction function. It should be given 
by 



(21) 



where dA^i/dv is the resonance profile approximation to 
dknijdv, i.e. 



dj^ 
dv 



|(n/||r||2p)r 



19683(2/ + l)7^a^ [y - vi^^^f 



[This equation is derived from Eq. (B5) of Ref. [1^, tak- 
ing only the leading-order (i'— J^Lya)"^ term and recalling 
that (2p||r||ls) = ^^^I'^a^ji^l'^ and t-Lya = There 
should technically be stimulated emission factors of I-I-/1. 
and 1 -|- /Lya in Eq. (|2T|). but in the vicinity of Lya 
< 10~^^ during our period of integration so we leave 
these out. Also £{i') should technically be computed in 
the rest frame of the hydrogen atom rather than in the 
comoving frame, but since Siv) varies extremely slowly 
with frequency (it does not possess a resonance at Lya) 
this correction is unimportant. 

Note that in Eq. (PT|) . the 2s level should be included 
in the numerator since it is possible for a 2s Is two- 
photon decay to produce emission within the frequency 
range of the Fokker-Planck code. In the blue wing of Lya 
> t^Lya), one should replace this with the 2s — > Is 
Raman scattering rate, 



dAo 



dv 



-(! + /.') 



dK2 
dv 



(23) 



[Note that there is no need to include 2s in the denom- 
inator, since the resonance approximation for its decay 
rate, Eq. ([H]), is zero.] Technically one should also in- 
clude continuum states in the sum, but since the true 
two-photon emission is dominated by decays from n > 3 
states rather than direct decays from the continuum we 
will not include the latter here. 



Equation (|2T|) is in general quite complicated. It can be 
evaluated under the approximation of Boltzmann equi- 
librium of the low-lying excited states. (The biggest ex- 
ception to this rule is the 2s : 2p ratio, which deviates 
from statistical equilibrium by 0.1% at z = 1190, 1% at 
z = 950, and 10% at z = 790. The 3s : 2p, 3p : 2s, and 
id : 2p ratios remain in Boltzmann equilibrium to < 1% 
throughout at all z > 700.) This allows us to construct a 
function £{v, T). This can be split into two contributions 



£{v,T) = £2s{'y,T) + Snyaiiy^T) 



(24) 



coming from the 2s and n > 3 levels respectively. Note 
that at V — vi^ya we must have £n>3 1 and £23 ~* 0. 
The functions £{v,T) is fit to < 0.3% accuracy over the 
range < 0.01 and T < 4700 K, where 1} = v/vLya - 1, 
by 



£„>3{v,T) = e 



-5.4iJ 



(25) 



and 



£2sii^, r) = 92.5e 



6.01J 



,/izyH„/T|^|3 



1 



|gtf/i^Ly„/T _ i| 1 + 0.321e-''''p-/r ■ 

(26) 

[The first fraction in this equation is physically motivated 
by the low-energy photon phase space factor oc the 
thermal stimulated emission (1 -I- /) or absorption (/) 



phase space density 1 / 



- 1 1 , and the Boltzmann 



(22) enhancement of 71 = 2 relative to n = 3 levels e'"'""/"^. 



The second fraction takes into account the fact that some 
of the Lya emission is preceded by Hp emission from the 
n = 4 levels, with 5"'"'^='°/^ representing the Boltzmann 
suppression of n = 4 relative to n = 3 hydrogen atoms.] 
The Voigt profile 4>\/{v) is computed using the integral 
formulation [38j . except in the far damping wings where 
we switch to the asymptotic expansion for \v — i^Lyal 3> 



2. Absorption 

True absorption is the inverse process of true emission, 
so the same matrix element applies to both cases. In 
particular, the ratio of two-photon absorption from Is to 
a given energy level u is related to the rate of emission 
via: 



N 



ab 



Qu'^lsfyfy' 



(27) 



where and f^' are the phase space densities associated 
with the two photons. We take v to represent the fre- 
quency of the photon near Lya and v' to represent that 
of the low-frequency photon. Since the lower-frequency 
photon comes from a blackbody distribution, we have 



1 + /. 



(28) 
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Further assuming that the u level is in Boltzmann equi- 
librium with 2p (a good approximation for the n < 4 
levels), we have 



92p 



iV|ab 3X1,/^, 



(29) 



(30) 



Using conservation of energy to find that hv' = £"„ 
-Els — we can simplify this to 



N\ 



ab 



X2p ^-h(v~v^^^-)IT, 



(31) 



This ratio applies for all of the excited states u, so it must 
apply to the total true emission and absorption rates as 
well. Thus we can solve for A^i|ab: 

iV.Ub = -^^e"(''-''--)/^^iV.|c„.. (32) 

We can simplify this further by defining the equilibrium 
number of photons per bin at line center, 



_ SttA In ly X2p 

*nn — 



(33) 



from which we convert Eq. ([5^ into 



N.U = f ^) e^^---- )/^'iV,|e„.. (34) 

Using Eq. (pO|) . we arrive at 



iVe, 



(35) 



The value of A^eq is provided by the interface f ^IIIBp . 



3. Scattering 

We now consider the change in frequency of photons 
due to resonant scattering off of hydrogen atoms, Eq. 
The typical fractional change in frequency is roughly 
Vth/c where Vth is the rms thermal velocity of the hy- 
drogen atoms. Nevertheless in a very optically thick line 
the net effect of many scatterings on the line profile may 
be important. We therefore write the Lyg transport in 
terms of a Fokker-Planck operator [H, [l^, [H, [4]| . 
In formulating such an operator, it is essential to be sure 
that the scattering term exactly conserves photons and 
preserves the equilibrium distribution oc e"'"'/-^™ [ioj . 
even after discretization. 

We define Fi to represent the net flux of photons from 
bin i -|- 1 to bin i due to scattering. In this way we have 



(36) 



This formulation guarantees exact conservation of pho- 
tons even in the discretized problem. In Fokker-Planck 
problems the fluxes are linear in the number of photons 
and its frequency derivative, i.e. in Ni and A^i+i — A^;, so 
we write 



F, 



(37) 



where Q and 77, are coefficients to be determined. In 
the equilibrium modified blackbody distribution, and for 
logarithmically spaced bins Aui cx Ui, we have 



(38) 



so in order for this to give zero net flux, we must have 

,3 



^1 



(39) 



Equation (|39p provides one constraint for two free pa- 
rameters Ci and r\i . The other constraint must come from 
fixing the diffusion coefhcient T>(y) to the correct value. 
We see that 

7V,|,e = C»-iA^»-i - (0 + m-l)N^ + (40) 

Taylor-expanding A'^; allows us to write 



1/ . ^d'^N, 



0th, 1st derivatives, (41) 



so the diffusion coefficient is ^{rji -f Ci-i) bin^ s ^. This 
can be written in the usual units of Hz^ s~^ by multiply- 
ing by the square of the bin width. 



(42) 



We then compare to the actual diffusion coefficient [4l| 

X>(l/) = i/l'LyaCr^TLya/s0v('^), (43) 

where a1 = v^^^^T^/ {niYic^) is the variance of the 
Doppler shift distribution due to motion of H atoms. The 
fraction of Lya absorptions that result simply in scatter- 
ing (as opposed to true absorptions) is /g; it is close to 
unity throughout the calculation, but a correct value is 
provided by the interface. 

Since and rji are slowly varying functions, we may 
replace Q-i — > Q in Eq. ([42|) and get 



■7^{m + Q = 



or 



?7. + 



i+1 



(44) 



(45) 



This and Eq. (|39p are sufficient to determine Q and rji. 

We handle the boundary conditions by disallowing any 
diffusion flux at either the red or blue boundary: F^i = 
Fm-i = 0, where M is the number of bins. 
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The line profile for scattering, (j)y{v), is in principle 
modified by the existence of neighboring resonances such 
as Ly/3. However, comparison of the Voigt profile to the 
actual cross section for Is — > Is scattering shows errors of 
< 2% in the frequency range of interest \'d\ < 0.01. Since 
the Is Is scattering makes only a < 0.45% correction 
to the recombination history, we ignore the "correction 
to the correction." 



4-. Hubble expansion and integration algorithm 

For logarithmically spaced bins in frequency, it is easy 
to compute the effect of the Hubble expansion: when 
the Universe expands by an amount A In a = A In j/, all 
photons simply shift into the next lowest frequency bin. 
Thus to account for the Hubble expansion, the contents 
of each bin are shifted down by one frequency bin at each 
time step: 

TVi =Ar,+i (previous). (46) 

Since there are only a finite number of frequency bins, the 
values of the number density of photons that redshift into 
the highest-frequency bin must be determined by some 
other means. This value is denoted A^in, and is provided 
by the interface ( mil B[) . This method of manually shift- 
ing photons to the left requires logarithmically spaced 
frequency bins and time steps. Moreover, the resolution 
Alnv is tied to the time step. (Despite this restriction, 
this method has the advantage of avoiding spurious nu- 
merical diffusion, which would arise if the Hubble expan- 
sion term were simply written as a differential operator 
with a discretized derivative.) 

Our method of solving Eq. psp is thus to apply an 
implicit ODE solver (backward Euler) to the emission, 
absorption, and scattering terms in Eq. p5p . evolve for- 
ward one time step, and shift the photons according to 
Eq. We repeat this basic operation until we reach 

the desired final redshift zumd- 

The abundance of photons in the highest frequency 
bin, Nn-i, is not specified by the above algorithm. Phys- 
ically it is determined by the phase space density of pho- 
tons redshifting into the line, in accordance with Eq. ([14]). 
This depends on the two-photon radiative transfer calcu- 
lation and is provided by the interface code. 

The emission, absorption, and scattering terms in the 
above equation can be written as a matrix equation: 

^i|cm+ab+sc = dij Nj , (47) 

where Cij is a tridiagonal matrix. We step forward using 
a backward Euler method: 

M±^1^.C.,(, + A,W« + A,); ,48) 

this is a first-order method but this is sufficient because 
of the extremely small time step. Inspection of the emis- 
sion, absorption, and scattering terms shows that Cij is 



tridiagonal and hence Eq. ([48]) is a tridiagonal linear sys- 
tem for {Ni{t + At)}. The M x M tridiagonal system 
can be solved by the usual 0{M) complexity method of 
using the i = equation to eliminate No{t + At), then 
using the i = 1 equation to eliminate A^i {t + At) , and so 
on. 



B. Interface 

The evolution equation for the {Ni} is only part of the 
recombination problem; it must interface to the multi- 
level atom code with the proper boundary conditions. 
This problem is considered here. The basic approach 
is iterative: the multi-level atom code is run first, to 
generate a table of input data for the Fokker-Planck code. 
Then the outputs of the Fokker-Planck code are used to 
apply corrections to the multi-level atom code, and so on 
until convergence is reached. 

The data passed from the multi-level atom code to the 
Fokker-Planck code at each time step are: 

1. The matter and radiation temperatures. 

2. The 2p state width (inverse lifetime including 
all processes that depopulate the state, including 
Lya decay, and bound-bound and bound-free ab- 
sorptions); 

3. The true emission rate of Lya photons H(a), com- 
puted using Eq. p^ . 

4. The equilibrium abundance of Lya photons 
Ncq/ Alnv, computed using Eq. ([55]) . 

5. The abundance of photons redshifting into the 
Fokker-Planck grid region per H nucleus per Hub- 
ble time, i.e. iVj\f_i/A Inz^ at vm-i- 

6. The fraction /inc of Lya absorptions that result in 
true absorption instead of scattering. This is deter- 
mined by the branching ratios for transitions out of 
the 2p level of hydrogen. Note that /s = 1 — /inc- 

7. The optical depth to scattering in Lya photons, 

Thyafs- 

In order to complete the iteration cycle, the Fokker- 
Planck code must return the corrections to Lya transport 
to the MLA code. This is done in several steps. First, 
we turn off the two-photon transitions in the MLA code 
involving frequencies between vq and vm^i- Then we cor- 
rect the usual equations for the Lya line with correction 
factors that cause it to produce the same outputs (net 
2p — > Is decay rate and photon phase space density at 
i/q) as the Fokker-Planck code. Explicitly, the standard 
equation for the Lya decay rate is 

_ 8ttH X2p Nin[{iyLya/l^M-l)a] 
X2p^ls,std[a) — -3 — — j , 

(49) 
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where the rate of incoming photons is measured at the 
frequency i^m-i > v-Lya at an earUer time since these are 
the photons that will reach Lya line center at scale factor 
a. [Compare to Eq. p2p .] The standard equation for the 
rate at which photons redshift out of Lya is 



fhyo 



std 



(a) 



X2p 

ix 



Is 



We replace these with the equations 



and 



/Lye- (a) = 6 (a) /Lya-. std (a)- 



(50) 



(51) 



(52) 



The correction factors $i(a) and ^2(0) are determined 
by the Fokker-Planck code as follows. The net decay 
rate is X2p-,is{a) = X^fio ^ ^i|cm+ab+sc. This equation 
is unstable as written if Ni is determined by plugging 
Ni into the evolution equations. The stable method is 
to find the change AiV in X]f=o ^ before and after the 
em+ab+sc time step, and write 



X2p^is{a) 



HAN 



(53) 



We can then find ^i(a) = X2p-,is{a) / X2p^is,stdia) . 

The red wing radiation correction factor ^2{o-) can be 
obtained by examining the phase space density of radia- 
tion emerging from the red wing of the line. In the ML A 
code with correction factor, this phase space density at 

will be 



6 (a 



X2pia) 
3xis{a) 



(54) 



The true phase space density is however known from the 
Fokker-Planck code, so one can solve for ^2(0)- 

Because ^i(a) and ^2(0.) are correction factors and are 
generally close to unity, we expect faster convergence by 
having the Fokker-Planck code return ^i(a) and £,2(0,) 
than absolute decay rates and phase space densities, so 
this is what we do. 

The frequency spacing and time step in the Fokker- 
Planck code are Alni/ = A In a = 8.5 x 10~^, which is 
5 times finer than the MLA code of Ref. [3 (A In a = 
4.25 X 10"^). Therefore the data provided by the MLA 
code are interpolated onto the finer grid required by the 
Fokker-Planck code. 

We find that only two iterations of alternately running 
the Fokker-Planck and MLA codes arc necessary. In our 
fiducial case, the first iteration leads to changes jAxel/a^e 
of at most 8.5 x lO"*^; the second iteration leads to a 
maximum change of 5.3 x 10~^; and the third iteration 
1.4 x 10-^. 

As a test, we have run the Fokker-Planck code with 
the scattering term Ni\sc turned off, and found agree- 
ment with the previous MLA code of Ref. [l^, with a 
maximum error \Axe/xe\ of 4 x 10~^ for 700 < z < 1600. 



Correction to recombination history from Lya scattering 
0.001 




-0.008 



700 800 900 100011001200 130014001500 1600 



FIG. 2: The correction to the recombination history due to 
LyQ diffusion. The numerical computation of i jllll is shown 
with a solid line, and i JIVI with a dashed line. Recombination 
is accelerated due to the additional redshifting of photons via 
atomic recoil. 



Results 



The results from the Fokker-Planck code are shown in 
Fig. [21 This run began at Zinit = 1605.5. As expected 
from heuristic arguments f ^II C[) . the rate of recombina- 
tion is accelerated by the inclusion of Lya diffusion. 

We have tested the convergence of our result with 
respect to the key numerical parameters. For exam- 
ple, if we only include scattering within ±500 bins of 
the line center instead of the full ±1000 bins, we find 
a maximum change in the ionization history jAi'e/xel 
of 10^^. As an additional test, we tried using a 2.5x 
coarser frequency binning for the diffusion code (so that 
the diffusion code takes 2 instead of 5 time steps in 
each step of the MLA code). The frequency range re- 
mained the same, so this corresponds to parameters 
A In J/ = A In a = 2.125 x lO'^ and M = 801 bins. This 
modification leads to a maximum change in the ioniza- 
tion history jAxg/xel of 5 x 10~^. 

The correction factors ^1(2:) and £,2(2) are shown in 
Fig. El 

We also show the radiation spectrum in Fig. |4l The 
radiation phase space density outside the diffusion code 
boundary (i.e. v < vq ov v > vm-i) is obtained from 
the MLA code with virtual levels as in Ref. [l^ , whereas 
for < v < I'M -I we have used the phase space density 
from the diffusion code. We show the no-diffusion case 
(MLA + virtual levels only) with the dashed line. The 
most noticeable effect of the diffusion is the increased 
intensity at v > ^'Lyc^ due to Lya photons diffusing to 
the blue side of the line. There is also an enhancement in 
the number of photons redshifting out of the line, which 
directly affects the recombination rate. 
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Correction factors, 5-| and 5; 




700 800 900 1 000 1 1 00 1 200 1 300 1 400 1 500 1 600 



FIG. 3; The correction factors ^i{z) and £,2{z)- The "stan- 
dard" result for infinitesimaUy narrow Lya hne with no damp- 
ing wing or diffusion effects is ^1(2;) = ^2(2) = 1- Note that 
over most of the recombination history these are greater than 
1, implying a faster 2p —r Is decay rate but also more photons 
redshifting out of the Lya line. The latter effect will lead to 
more two-photon absorption at low redshifts. [The glitch in 
^i{z) at 2; « 1360 is a startup transient from Ly/3 at z = zinit 
redshifting into Lya; given that the overall effect of the scat- 
tering correction is < 0.6% in d, this is far too small to affect 
CMB results.] 



Lya line profile at z=1 006 

^(A) 

1240 1230 1220 1210 1200 




I ' ^ ' — ' 1 

0.735 0.74 0.745 0.75 0.755 0.76 0.765 

V (Rydbergs) 



FIG. 4: The radiation spectrum in the vicinity of the Lya 
hne at z — 1006. The solid line shows the results of the 
diffusion code, and the dashed line shows the old code with no 
diffusion [l^. The vertical dotted lines show the boundaries 
of the Lya diffusion region. The "chemical equilibrium" curve 
shows the approximation ~ (a;2p/32;is)e~'^*'^~'^'"y°-'^"^' that 
should be valid in the immediate vicinity of line center. Note 
that i/Lya ~ 0.75 Rydbergs. 



escape rate. This approach is a valuable complement 
to the fully numerical method: it contains additional 
approximations, but it provides a better understanding 
of the physics and a check of the much more sophisto- 
cated Fokker-Planck code/interface. Analytic corrections 
(or more precisely, corrections based on reduction of the 
problem to a simple ODE) have been used by previous 
authors [34|. j35, -12] . The specific implementation here is 
an extension of the two-photon analysis of Ref. [l^ to 
include Lya diffusion. 

In Ref. the problem of emission and absorption in 
the Lya damping wings was considered with no diffusion 
(i.e. without considering the change in frequency during 
a Is ^ Is scattering). The red and blue damping wings 
were handled separately since the radiative transfer phe- 
nomenology is different. In both cases, the equation for 
the radiation field is written down and is approximated 
by its time-steady form. Then: 

1. In the red wing, wc find the correction to f^, 
and based on the concept of the fiux of photons 
[Eq. (|lip ] we find a correction to the rate of Lya 
escape. Since the correction is always positive 
[Eq. ([T^ ] the red wing corrections always increase 
the recombination rate, leading to lower ionization 
fraction. 

2. In the blue wing, we find the number x+{t) of spec- 
tral distortion photons in the blue wing of Lya per 
H nucleus according to the time-steady equation. 
This function starts at zero before recombination, 
reaches a positive maximum, and then declines to 
zero at late times as all photons rcdshift to lower 
frequencies. An additional downward decay rate 
i2p^is = i+{t) is grafted on to the MLA code to 
account for the 2p — > is decays that are required 
early during recombination to build the spectral 
distortion, and then the excitations that occur later 
during recombination as this distortion redshifts 
into Lya. Note that this process accelerates re- 
combination at early times (i+ > 0) but delays it 
later {x+ < 0). 

We now extend the treatment of Ref. [19| in the red 
f ijIV Ap and blue (fJVB]) wings. In all cases we ne- 
glect the variation in phase space density (i.e. factors 
of v/vhya) across the Lya line. 



A. Red wing 



IV. ANALYTIC APPROXIMATION 



Without frequency diffusion, the radiative transfer 
equation, Eqs. (82, 85, 87) of Ref. [l|, is 



We now consider a completely different approach to the 
Lya diffusion problem, in which simple analytical calcu- 
lations are used to estimate the correction to the Lya 



(55) 
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where 



W 



47r2 



E 

nZ,n>3 



2Z + 1 



i,2p 



(56) 



is the width over which Lya is optically thick to true 
absorption. 

We can add the frequency diffusion to this equation by 
recalling the diffusion term [4l| , 



fl 

Of 



(57) 



[Note that the (h/T^fi, term accounts for recoil.] In the 
far damping wings, we have 



(58) 



(see Ref. |41| and use the damping wing approximation to 
the Voigt profile for large ly — Vhya ) ■ Adding this equation 
to Eq. gives 



Hv 



W 



dU 

diy {v ~ y-Lya) 

Cr^TLya/s^LyQ 



Xi'-i^Lyc)/T, 



X2p 



47r2 



d_ 

dv 



1 



As in Ref. [l^, we make the change of variables 



I^Lya 



-y 



and 



/. = ^$(y). 



(59) 



(60) 



(61) 



This, combined with dropping the fi, term (time-steady 
approximation) and taking « Tj- (appropriate dur- 
ing the recombination era for the purposes of computing 
small corrections) simplifies Eq. ([59|) to 







dy 



W 

-^(e^$ 



S- 



y 



dy 

where W = hW/T, and 



dy 



47r2T3 



, (62) 



(63) 



This results in a dimensionless equation, Eq. (|62p . which 
depends on two constants W and 5'. The constant W 
determines the strength of the true absorption: the Lya 
line is optically thick to true absorption out to frequencies 
vi^yaiWT^/h. The constant S quantifies the importance 
of frequency diffusion relative to Hubble redshifting at 
frequencies i^Lya ± Ti-/h. In practice both are ^ 1. 



(a) Correction to Lya escape rate 



0.1 



CD 



0.01 




0.001 



(b) Dimensionless abundance of photons in blue wing 



0.1 



0.01 > 



0.001 




FIG. 5: The analytic correction factors x(^) (top panel) 
and J{W, S) (bottom panel) associated with transport in the 
Lya line. The thick bold line shows the factors without fre- 
quency diffusion, i.e. for S = 0. The thin lines show the 
factors for S = 10"^ 10"^ 10"^ 10"", and 10"^ from bot- 
tom to top. Note that both emission/absorption in the far 
damping wings (parameterized by W) and frequency diffu- 
sion (parameterized by S) tend to increase the escape rate 
and the number of distortion photons in the Lya blue wing, 
as expected. 



Our time-steady equation, Eq. (j62p . is very similar to 
Eq. (93) of Ref. [19j , and it satisfies the same boundary 
condition: $ = 1 at y = 0, since at line center we reach 
equlibrium and have /i/Ly„ = a^2p/(3xis). The second- 
order differential operator complicates the solution and 
necessitates an additional boundary condition that <f> not 
diverge as i/ ^ — oo. The numerical solution is presented 
in Appendix [XI Just as in Ref. [l^, the correction to the 
net 2p — > Is decay rate is 



Ail 



^Lya 
Thya 



(64) 



where x = ^{y = —oo). The correction x — 1 is now a 
function of the two dimensionless constants, W and S. 
It is shown graphically in Fig. [51 
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B. Blue wing 

The frequency diffusion in the blue damping wing leads 
to a modification of the number of spectral distortion 
photons x^{t) per H atom in the blue wing of Lya. As 
in Ref. [l^, this can be approximated as 



x+ 



(65) 



where y is dimensionless frequency and A/^ is the dis- 
tortion contribution to the phase space density. Writing 
the spectral distortion as 



(66) 



Rcf. [19| showed that in the absence of frequency dif- 
fusion the rescaled spectral distortion ^{y) satisfied the 
equation 



dy y^ 



(67) 



in the time-steady approximation, i.e. the same equa- 
tion as occurs in the red wing. [This is Eq. (109) in 
Ref. [T5|; the missing y^ in that paper is a typo.] The 
inclusion of frequency diffusion proceeds exactly analo- 
gously to i]IV Al yielding 



= 



dy 



dy 



y 



dy 



(68) 



The abundance of distortion photons in the blue wing is 
then 



where 



^^'^lygT, / X2p 



IiW,S) 



-'^-^y^/TAj(^W,S), (69) 



'^iy)dy 



(70) 



is a dimensionless integral. Values of T{W, S) are com- 
puted according to the method in Appendix [X] and plot- 
ted in Fig. El 



C. Implementation and results 

Following the approach of Ref. [l^, we first compute 
the Lya transport parameters W{z) and S{z) for the 
pure MLA code with all two-photon transitions and scat- 
tering effects turned off. We then turn on the analytic 
corrections in Ref. [l^ associated with the stimulated 
2s —^ Is decays and nonthermal absorption, two-photon 
decays from n > 3 levels, and Raman scattering. The 
two-photon decays from n > 3 levels depended on the 
function x{^) = ^{u = — oo|W^) for the sub-Lya decays 



Lya line parameters 



CO 

b 




700 800 900 1000 1100 1200 1300 1400 1500 1600 
z 

FIG. 6: The dimensionless Lya transport parameters W{z) 
and S{z). 



(i.e. those in which both of the emitted photons have 
< i^yai which usually means one photon emerges in 
the red damping wing of Lya) and I{W) for the super- 
Lya decays (where one photon has v > Vhya, usually 
in the blue damping wing of Lya). We can account for 
scattering semianalytically by replacing x(VF) and T{W) 
with their generalized values xi^j 5*) and T{W, S) de- 
rived here. We may then compare the resulting recombi- 
nation histories with and without Lya scattering. 

The transport parameters W{z) and S{z) are shown 
in Fig. E 

The correction to the recombination history can be 
obtained by comparing the "old" Xe{z) using x(^) 
and 1{W) without scattering to the "new" Xe{z) using 
x{W,S) and I{W,S). The correction is shown by the 
dashed line in Fig. [2] Note the qualitative agreement 
with the fully numerical result, although at low redshifts 
our analytic approximation overestimates the correction. 



D. Comparison with other computations 

The changes in the recombination history that we have 
found amount to corrections of at most 0.45%. This is 
less than the correction computed by several other au- 
thors. This section discusses some possible explanations 
for the apparent discrepancies. In some cases, the expla- 
nation lies with the fact that the corrections to the re- 
combination history (or the effective escape probability) 
from emission/absorption in the damping wings are not 
additive with those from scattering - the net effect of in- 
cluding both is less than one would expect from adding 
the contributions to Axe/xg. Also previous results on 
Lya transfer did not include the deviation of emission 
versus absorption profiles, which can have a major im- 
pact on the results - e.g. without this we would have 
X{W, S* = 0) = 1 for any W^. 

Chluba & Sunyaev |35| consider the time dependence 
of the radiation intensity in the Lya line (i.e. non-quasi- 
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stationarity) and the consequent effect on recombina- 
tion. They include true emission and absorption in the 
damping wings but not scattering, and so in terms of 
the physics their result is most comparable to the treat- 
ment of two-photon decays by Hirata p^ . In partic- 
ular, Hirata found that the dominant time-dependent 
correction was that associated with the blue damping 
wing of Lya, i.e. with the time dependence of i+ (dis- 
cussed here in ijIVB[) . Both Chluba & Sunyaev [S^ (see 
their Fig. 12) and Hirata [l^ (see his Fig. 8) find that 
the non-quasi-stationarity leads first to an accelerated 
recombination and then a delayed recombination as the 
spectral distortion redshifts through Lya, but Chluba & 
Sunyaev find an effect up to a factor of ^ 3 larger. We 
suspect this is due to their neglect of the deviation of 
emission versus absorption profiles, i.e. the factor in 
Eq. dMl)- Without this factor (and with S* = 0) we de- 
rive the solution ^'(j/) = 1 — e~^/^ and hence the integral 
T{W, 5 = 0) = '^Jj^dy = oo, and so the analytic ap- 
proximation in Ref. [l9| would yield an infinite correction 
in this approximation. Chluba & Sunyaev find a finite 
correction because of the finite duration of recombination 
(they use an exact treatment of non-quasi-stationarity 
rather than Hirata who treats the time dependence as a 
perturbation), but the effect is still large. 

Grachev & Dubrovich [sj computed a correction to 
Lya escape based on the modified escape probability of 
Grachev [44| . The latter did not include non-time-steady 
effects and was based on the rate of redshifting of Lya 
photons out of the line, similar to our ijIV Al They also 
did not include the deviation of emission versus absorp- 
tion profiles, which is equivalent to ignoring the in our 
Eq. (p^ . This led them to the analytic approximation of 
Grachev [431, which is 



Change in power spectra 



X - 1 = PG 



1 



= (4- 
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(71) 



3(2 + 4) 

where Grachev's dimensionlcss parameters are CTq = 
32/3^^/^ and pG = [Grachev [H uses i{-oo) 

in place of our x-] We have integrated our equation by 
the method of Appendix |^ without the factor and 
find that this agrees with Eq. ((7T|) to better than 10% 
in the range of interest. We note that our change in x 
due to scattering, xi^^S) — x(W, 0), is increased if we 
turn off this term, in qualitative agreement with the fact 
that Grachev & Dubrovich [IJl find a larger change in 
the recombination history due to Lya scattering. 



V. IMPLICATIONS FOR CMB ANISOTROPIES 

In Fig. [71 we evaluate the effect of the Lya diffusion 
correction on the CMB anisotropics, i.e. 

ACj^ Cj^(with scattering) 



C 



TT 



(without scattering) 



1. 



(72) 



Note that the two-photon transitions are turned on in 
both the "with diffusion" and "no diffusion" case. The 
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FIG. 7: The change in CMB power spectra for the temper- 
ature (top panel), polarization (middle panel), and the cor- 
relation coefficient (bottom panel) due to inclusion of Lya 
scattering. The overall tilt and oscillating behavior are both 
the result of a faster recombination and hence higher redshift 
of the surface of last scattering: the tilt from the reduced 
Silk damping, and the oscillations due to the smaller acoustic 
horizon. 



CMB power spectrum is computed using the Boltzmann 
code Cmbfast [1]. The most obvious effect is the oscil- 
lation in IS.Ci/Cf. these are due to a slight shift in the 
acoustic scale to higher £ since the faster recombination 
results in a higher-redshift surface of last scattering. The 
overall tilt is partly a result of the reduced Silk damping: 
the faster recombination gives less time for the acoustic 
oscillations to be damped by photon diffusion, and hence 
the small-scale perturbations are not as suppressed as in 
the standard scenario. Also the reduced electron density 
at 2: ~ 900 implies a lower optical depth after the surface 
of last scattering and hence less washing out of the small- 
scale features in the CMB. This latter effect is slightly 
overestimated by the analytic computation, which is why 
the analytic result (dashed line in Fig. [71) has a slightly 
larger high-£ power spectrum. 

As in Ref. [l^, we may evaluate the importance of 
the Lya scattering correction for a given experiment by 
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TABLE I: The magnitude of the correction Z (in number of 
sigmas) to the CMB power spectrum for several CMB exper- 
iments. This is shown both for the Lya diffusion correction 
(center column) and for the combined two-photon decay cor- 
rections [l9| and diffusion correction (right column). 

Experiment Z Z 

scattering 27-(-scattering 

WMAP 5yr 0.07 0.23 

ACBAR 2008 0.04 0.17 

Flanck 0.92 5.01 



ing from zero to 0.4% at £ < 3000, of the order of O.Qcr 
for Planck. 

Our ultimate goal in hydrogen recombination is to de- 
velop a complete theory with an error budget as has been 
done for helium [2^. This will be based in part on the 
work presented in this paper, but will also require investi- 
gation of many minor radiative and coUisional processes. 
Eventually the theory must be encapsulated in a code 
fast enough to use in Markov chain parameter estimation 
techniques, which could be the result of either analytic 
simplifications or interpolation codes . 



considering 



Z' 



(73) 



here Fui is the experiment's Fisher matrix and Z is the 
maximum number of sigmas by which any parameter fit 
could be affected. (An individual cosmological parameter 
may or may not be affected, depending on whether its ef- 
fect on the CMB power spectrum is similar to that caused 
by the inclusion of scattering.) We consider the impli- 
cations for the WMAP 2008 (5-year) data 3 (including 
beam and point source errors), the Arcminute Cosmology 
Bolometer Array Receiver (ACBAR) 2008 power spec- 
trum [ist (including beam and calibration errors), and 
the upcoming Planck data as forecast using the noise 
curves for the 70 GHz (Low-Frequency Instrument) and 
100 and 143 GHz (High- Frequency Instrument) channels 
in the Blue Book [4j and assuming a usable sky coverage 
of /sky = 0.7. 

Results of this analysis are displayed in Table HI It is 
seen that the correction due to Lya scattering is small - it 
is only 0.9(7 for Planck. Moreover, it goes in the opposite 
direction to the two-photon corrections found by Ref. [ll| 
(the scattering raises the high-^ power spectrum whereas 
the two-photon corrections lower it). 



VI. DISCUSSION 

We have examined the effect of multiple scattering 
on the H i Lya escape problem during the cosmolog- 
ical recombination epoch. Both numerical and analytic 
tools were developed to treat this problem simultaneously 
with emission and absorption in the damping wings. We 
find that scattering increases the escape probability and 
speeds up recombination because atomic recoil leads to 
a systematic shift of the photons to the red wing of the 
Lya line. While this is in qualitative agreement with 
previous results [H, [13] , the magnitude of the scattering 
correction is less than previously suggested. We believe 
this is because the proper treatment of emission and ab- 
sorption leads to a smaller effect from scattering. The 
modified treatment of the H i Lya resonance results in a 
very small correction to the CMB power spectrum, rang- 
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APPENDIX A: SOLUTION TO TIME-STEADY 
DIFFUSION EQUATION 

In this appendix, we consider the solution to the 
dimensionless time-steady radiative transfer equation, 
Eq. (|62p . In the red damping wing, y < 0, we desire the 
value of the solution at large negative values <&(— oo). In 
the blue wing, y > 0, we desire the integral <^{y)dy. 
Our ODE is different from the case of no frequency diffu- 
sion [l^ because the diffusion operator renders it second- 
order. It can be solved by defining the variable 



(Al) 



Here S(y) can be thought of as a dimensionless flux of 
photons passing through y due to frequency diffusion 
(Sy^^d^/dy) and recoil (Sy^^^). It is exactly zero with- 
out the diffusion/recoil terms. This leads to the linear 
system 



dy 
dE 
dy 



— S — <I> 

s ' 

Wy~^{ey<^-1)-^E 



(A2) 



which is singular at j/ = 0. 

We then implement a shooting method in the red 
damping wing, and a numerically stable modification of 
the shooting method in the blue wing. 
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1. Red wing 

First consider the red wing, y < 0. The boundary 
conditions are that $(y = 0) = 1 and ^{y ~ — oo) is 
finite. In order for = — oo) to be finite, inspection 
of the first equation in Eq. (jA2p shows that we must 
have S(j/ = — oo) = 0. Therefore given any choice of 
X = ^{y ~ ~oo), wc may construct a solution $(y|x) to 
Eq. (|A2p by taking $ = x and S = at large negative 
y and integrating toward y = using a second-order 
implicit ODE integrator. 

The problem remains to choose the correct value of 
X- In general for small negative y, the behavior of the 
solution is that <f>(y) remains finite because of the y^ 
suppression in the first equation of Eq. (|A2[) . while 



can have a 



y 



divergence as one can see from the 



second equation: 



lim $(?/) 



W{1 - a) 



y 



0{y°). (A3) 



The desired solution is that corresponding to the bound- 
ary condition a = 1. As we integrate the solution $(y|x) 
toward y = 0, we can determine the value a as a function 
of X- Since the ODE is linear, a(x) is a linear function 
of X, and it suffices to obtain a(x = 0) and a{x = !)■ 
Then the desired value of x that satisfies the line center 
boundary condition can be obtained by linear extrapola- 
tion: 



X{a = 1) = 



1 - a(0) 
a(l) -a(0)' 



(A4) 



Since the fundamental problem is to find $(y = — oo) 
we may simply report x(a = !)■ 



2. Blue wing 

The solution in the blue wing, y > 0, must satisfy the 
boundary conditions <&(?/ = 0) = 1 and = +oo) = 0. 
As before, we will implement a shooting method in the 
+y direction, i.e. starting at small positive y ~ e and 
integrating toward +oo. But here the numerical calcula- 
tion is much trickier: whereas in the red wing the grow- 
ing mode of Eq. (|A2p had only a power-law divergence 
(S oc y^^) as one approached line center, here the grow- 
ing mode grows faster than exponentially. Therefore if we 



push to large positive y an overflow error occurs. More 
seriously, this growing mode implies that the correct ini- 
tial condition E{y = e) cannot be well represented even 
in double precision. 

A slow but effective way to solve this problem is to 
write a function that takes in an initial value of yi and 
a dimensionlcss photon density "^(j/i) at that point, and 
find the critical value Sc[yi, ^•(i/i)] of S that satisfies the 
boundary condition at y ~ +oo (i.e. is non-divergent). 
Here Sc can be determined by a shooting method using 
a second-order implicit ODE integrator with two trials, 
taking advantage of linearity, just as we did for the red 
wing. We also determine how far one can integrate before 
the dangerous growing mode kicks in. This can be found 
from the function 



F[?/|yi,$(yi) 



a$[y|?/i,$(?/i),^(?/i)] 



9S(yO 



oV-Vi 



(A5) 

which takes the initial value F{yi) = and rises toward 
oo as 2/ increases. Wc solve for i/thbi, ^(j/i)] at which 
F = 10^, i.e. where the growing mode has incresed by 4 
orders of magnitude. The exponential e^~^' is included 
because the physical solution has ^ dependence at 
large y. 

We then begin our solution to Eq. (jA2p by starting at 
y — t, ^{y) — e~'^ , and S(y) = Sc(e, e^^). The implicit 
ODE integrator is used to integrate to larger y, but only 
until we reach ?/th- Beyond this point the growing mode 
is dangerous, so we recompute the critical Sc at [y, <&(?/)], 
take this as a new initial condition, and keep integrating 
until the new i/th- This procedure "resets" the unstable 
mode each time a threshold value ?/th is reached. 

At very large values of j/, where the growing mode 
grows very quickly, the above procedure becomes very 
slow. Fortunately, we find that at this point $ is well- 
represented by the large-y asymptotic result <&(?;) 
e~^. We may then extract the integral XiW.S) — 
/„°°$(y|M/,5)d2/. 

In the real Universe, the solution at very large values 
of y is not accurate because of the presence of other res- 
onances (e.g. Ly/3 is at ?/ ^ 8 at z = 1000). Fortunately, 
these large values of y do not contribute significantly to 
the integral T{W^ S) because of the suppression of 
the integrand; it is for this reason that the analytic ap- 
proximation based on I(W, S) works so well. 
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